The Geometrical Structure of Disordered Sphere Packings 
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The three dimensional structure of large packings of monosized spheres with volume fractions 
ranging between 0.58 and 0.64 has been studied with X-ray Computed Tomography. We search 
for signatures of organization, we classify local arrangements and we explore the effects of local 
geometrical constrains on the global packing. This study is the largest and the most accurate em- 
pirical analysis of disordered packings at the grain-scale to date with over 140,000 sphere coordinates 
mapped. We discuss topological and geometrical ways to characterize and classify these systems, 
and discuss implications that local geometry can have on the mechanisms of formation of these 
amorphous structures. 

PACS numbers: 45.70.-n Granular Systems 45.70.Cc Static sandpiles; Granular Compaction 45.70.Qj Pattern 
formation 
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I. INTRODUCTION 

When balls are poured into a container they arrange 
themselves in a disorderly fashion with no obvious sym- 
metries or repetitive patterns. However, disorder does 
not mean randomness. Indeed such systems are lo- 
cally highly structured in a hierarchical organization 
which tries to achieve the goal of maximal compaction 
under the unavoidable geometrical constraints of non- 
interpenetration, satisfying simultaneously the condition 
of force and torque balance on each ball. This leads to 
very complex structures which show signs of organiza- 
tion but nevertheless have so far eluded all efforts for 
a simple and clear classification. In order to fully clas- 
sify the state of a disordered system, such as a granular 
packing at rest, the exhaustive details about the exact 
position, orientation and shape of each grain is, in prin- 
ciple, needed. However, part of such information is at 
best redundant or even irrelevant and several degenerate 
states with different microscopic realizations can share 
the same macroscopic properties. To determine which are 
the accessible configurations at the local level, and to un- 
derstand which are the possible combinations which gen- 
erate the global packing is of fundamental importance. 
Indeed, finding measures for the local and global, hierar- 
chical organization is the essential starting point towards 
the understanding of the basic mechanisms which form 
these structures. It is also an essential step in the de- 
velopment of technologies which enable us to control and 
'tune' the structure of amorphous materials. 

Until now the empirical investigation of the geometri- 
cal structure of these systems have been limited by the 
very little availability of accurate experimental data. In- 
deed, after the seminal works of Bernal, Mason and Scott 
HI S 01 > it nas been only very recently that the use of 
tomography has allowed to investigate three dimensional 
structure from the grain level up to the whole packing. 

The first work which uses tomographic techniques de- 
voted to the investigation of granular packing is by Sei- 
dler et al. in 2000 Other two works by Sederman 



et al. [j| and Richard et al. [g followed respectively 
in 2001 and 2003. Confocal microscopy techniques have 
been also used to reconstruct 3D images of a dense pack- 
ing emulsion of oil droplets and to count contacts in 
glass beads However, all these works concern rather 
small sample sizes and focus the analysis on few partic- 
ular topics. In this paper we present an empirical inves- 
tigation by means of X-ray Computed Tomography on 
very large disordered packings of monosized spheres with 
packing densities ranging from 0.58 to 0.64. (The pack- 
ing density is the fraction of volume occupied by the balls 
divided by the total volume, and it is often called 'vol- 
ume fraction' in the literature). This study is the largest 
and the most accurate empirical analysis of disordered 
packings at the grain-scale ever performed. A packing 
realization is shown in Fig^ Preliminary results were 
presented in ref. [jj. Here we perform a more extensive 
and complete investigation using an improved algorithm 
to calculate the positions of the spheres. Additional ma- 
terial including the sphere coordinations of some of the 
samples can be found in [lOj. 

In order to help the reader this paper is organized in 
several sections each addressing different aspects. Each 
section has been designed to be as self-containing as pos- 
sible. Cross references among sections guide the reader 
who may wish to focus on specific topics. 

In Section [H] the experimental apparatus and the rel- 
evant methodology to extract geometrical information 
from the tomography data are described. 

In Section ITU the number of neighbors for each sphere 
in the packing is studied and a new tool to deconvolute 
the contribution of touching neighbors from the contri- 
bution from near neighbors is introduced. Implications 
on mechanical equilibrium are also discussed. 

In Section HVI the hierarchical structure of the contact 
network is analyzed in terms of a shell map [Tl"j |. 

The local symmetries are explored in Section by us- 
ing spherical harmonics decomposition [T^ . 

In Section IVII the compactness of the average local 
packing is discussed and compared with crystalline pack- 
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FIG. 1: A reconstruction of a packing of ~ 150, 000 spheres 
in a cylindrical container (sample C). 

ings. 

We analyze and discuss the results for the radial dis- 
tribution function in Section [VIII . 

The density fluctuations at sample level and at grain 
level are investigated in Section IVIIII The implication 
that sample geometry can have on the dynamical forma- 
tion of these systems is discussed in Section IIXI 

A conclusion summarizes the main results and perspec- 
tives. 



II. EXPERIMENTAL APPARATUS AND 
METHODOLOGY 

The empirical studies reported in this paper concern 
the analysis of 6 samples made of mono-sized acrylic 
beads packed in a cylindrical container with an inner di- 
ameter of 55 mm and filled to a height of ~ 75 mm. In 
particular we have: 

• Two large samples containing ~ 150, 000 beads 
with diameters d = 1.000 mm and polydispersity 
within 0.05 mm; 

• Four smaller samples containing ~ 35,000 beads 
with diameters d = 1.59 mm and polydispersity 
within 0.05 mm. 

An independent estimation of the polydispersity for 
the small spheres was performed by weighing 200 beads 
and computing their standard deviation. The estimated 
value for the relative statistical error on the sphere di- 
ameters is 1.5%. 



A. Sample Preparation 

The 6 samples (named hereafter with labels A-F) have 
been prepared at different packing densities p ranging 



between 0.586 (sample A) to 0.640 (sample F). Tabled 
reports in the second column all the sample densities. 
The two packings at lower densities (A, B respectively 
with p ~ 0.586 — 0.596) were obtained by placing a stick 
in the middle of the container before pourin g th e beads 
into it and then slowly removing the stick |l3j| . Sam- 
ple C (p = 0.619) was obtained by gently and slowly 
pouring the spheres into the container. Whereas, the 
sample D (p = 0.626) was obtained by a faster pour- 
ing. Higher densities, up to p ~ 0.63, were achieved by 
gently tapping the container walls. The densest sample 
at p — 0.64 was obtained by a combined action of gen- 
tly tapping and compression from above (with the upper 
surface left unconfined at the end of the preparation). 
To reduce boundary effects, the inside of the cylinder 
has been roughened by randomly gluing spheres to the 
internal surfaces. 



B. XCT Imaging 

A X-ray Computed Tomography apparatus (see Sakel- 
lariou et al. |14|) is used to measure the density maps 
of the samples. The two large samples (A, C) were anal- 
ysed by acquiring data sets of 2000 3 voxels with a spatial 
resolution of 0.03 mm; whereas the four smaller sam- 
ples (B, D, E, F) were analysed by acquiring data sets 
of 1000 3 voxels with a spatial resolution of 0.06 mm. 
After segmentation (see Sheppard et al the sample 
data sets are reduced to three-dimensional binary images, 
representing two distinct phases, one associated with the 
spheres and the other with air space. The effective spa- 
tial resolution of this technique is limited by the finite size 
of the X-ray source, surface scattering of the low energy 
X-ray and intrinsic blurring from reconstruction. From a 
careful analysis of the reconstructed samples we observed 
that the combination of all these factors generates some 
fuzziness in a region between one and two voxels around 
the sphere surfaces. 



C. Sphere Centres 

In order to proceed with the analysis of the geomet- 
rical and statistical properties, the position of all sphere 
centres are calculated from three dimensional binary im- 
ages. Our approach is to find the sphere centres by 
moving a reference sphere (S) throughout the packing 
(P) and measuring the local overlap between S and P. 
This corresponds to a 3-dimensional convolution: P * S. 
The regions with larger overlaps are the ones around the 
centers of the spheres in the packing. In order to iso- 
late these regions we apply a threshold on the intensity 
map resulting from the convolution. The centres of mass 
of such regions are good estimations for the positions 
of the packed spheres. This method is made numeri- 
cally highly efficient by applying the convolution theorem 
which allows to transform the convolution into a prod- 
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FIG. 2: Number of detected spheres in an internal region 
of samples A (top) and D (bottom) reported as function of 
the reference sphere radius and the threshold. Threshold is 
expressed in relative units with respect to the value of the 
highest intensity peak in that region. 



tion of the reference sphere size and the threshold. When 
the reference sphere is too small only one cluster is de- 
tected independently on the threshold. Similarly when 
the threshold is too small, peaks cannot be isolated and 
one spanning cluster is also observed. On the other hand 
there is a rather large region of the two parameters where 
the same amount of spheres is detected. We know that 
precision increases with the size of the clusters on which 
the center of mass is computed. Therefore the best choice 
of parameters is the one which leads to the largest clus- 
ters. This requires the smallest possible threshold and 
reference sphere compatibly with the correct detection 
of each sphere in the system. We chose the threshold at 
0.95 (samples A, C, E ) and at 0.93 (samples B, D, F), 
and we fixed the reference sphere radius at 13 (sample A, 
C), 11 (sample E) and 10 (samples B, D, F). Obtaining 
typical cluster sizes of ~ 80 voxels (samples B, D, E, F) 
and ~ 400 (samples A, C). This implies precisions on the 
sphere centers respectively within 3% and 0.5 % of the 
voxel sizes 561. 



D. Central Region 

In order to reduce boundary effects, all the analysis 
reported hereafter have been performed over a central 
region (G) at 4 sphere-diameters away from the sam- 
ple boundaries. Note that spheres outside G are consid- 
ered when computing the neighbouring environment of 
spheres in G. In Table [I] the number of spheres in this 
region (Ng) is reported for each sample. 



uct in Fourier space: T[P * S] = ^[P^IS], where F 
represents the (fast)Fourier Transform. The algorithm 
proceeds in 4 steps: 1) fast Fourier transform of the bi- 
nary image [T\P\)\ 2) transform the digitised map of the 
reference sphere (^[5]); 3) perform the direct product 
between these two; 4) inverse-transform of the product: 
^ r_1 [^ r [P]J r [S']] = P * S. The result is an intensity map 
of the overlap between the reference sphere and the bead 
pack, where the voxels closer to the sphere centres have 
a highest intensity. A threshold on the intensity map 
locates the groups of voxels surrounding the sphere cen- 
ters which become isolated clusters. The sphere centers 
are calculated as the centre of mass (intensity) of these 
clusters. 

The precision on the estimation of the sphere-center 
positions can be evaluated considering that the spatial 
resolution is within one and two voxels. Therefore the 
precision on the centre of mass of a cluster of v vox- 
els must be within 2/u. In the procedure to locate the 
sphere-centers one has two dajusTable parameters: the 
reference sphere size and the threshold. We have searched 
for the optimal choice of the parameters by varying these 
two quantities and computing the resulting number of 
spheres detected in a given portion of the sample. In 
Fig|21the number of detected spheres is reported in func- 



III. NUMBER OF NEIGHBORS IN CONTACT 
A. Near Neighbors 

Let us start the analysis of the local packing config- 
urations by exploring local neighborhoods. The average 
number of spheres in contact with any given sphere is the 
primary and most investigated parameter in the litera- 
ture on granular packings H 0, EL IE H EH El 
0,0. Indeed, this is a very simple topological quantity 
which gives information about several important proper- 
ties of the system. Unfortunately, although simple in its 
definition, such a number is an ill-defined quantity from 
a experimental point of view. The reason is that, from 
a geometrical perspective, the information about the po- 
sitions, and eventually the sizes of all spheres is not suf- 
ficient to determine such a number: two spheres can be 
infinitesimally close but not in touch. In the literature 
several physical methods have been used Q, H, 0] , but 
they all encounter problems essentially associated with 
the uncertainty in the threshold distance used to define 
the maximum allowed gap between apparently 'touching' 
spheres. 

An exact computation of the number of touching 
spheres from a geometry alone is, in general, an impos- 
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density N Ng n c 


nt(l) 


m(1.02) 


m(1.05) 


m(l.l) 


t 


A 0.586 ± 0.005 102,897 54,719 5.81 


3.0 


5.5 


6.7 


7.5 


0.014 


B 0.596 ±0.006 34,016 15,013 5.91 


2.9 


5.9 


6.8 


7.7 


0.011 


C 0.619 ±0.005 142,919 91,984 6.77 


3.5 


6.4 


7.5 


8.4 


0.013 


DO. 626 ±0.008 35,511 15,725 6.78 


3.3 


6.0 


7.5 


8.4 


0.017 


E 0.630 ±0.01 35,881 15,852 6.95 


3.4 


6.3 


7.6 


8.6 


0.016 


F 0.640 ±0.005 36,461 16,247 6.97 


3.3 


6.9 


7.9 


8.9 


0.011 



TABLE I: Sample density and their interval of variations (±) within each sample; number of spheres in the sample (N); number 
of spheres in the central region (Ng); estimated average number of neighbors in contact (n c ), average number of neighbors 
at given radial distance (n t (r) with r = 1,1.02,1.05,1.1 diameters). Standard deviation (£) calculated from the probability 
distribution for radial distances smaller than d between pair of sphere centers. 




FIG. 3: (symbols) Behavior for the average number of sphere 
centers within a radial distance r. (lines) Complementary 
error function n^r)* 1 ' (Eq0 normalized by best-fitting the 
agreement with the data in the region r < d. The averages 
(d) and the standard deviations (£) are calculated from the 
probability distribution for radial distances smaller than d 
between pair of sphere centers. The re-normalized comple- 
mentary error function fits well the data for r < d. After 
this value near-neighbors not in contact start to contribute 
significantly to n t (r) and the two behaviors split. The 'de- 
convoluted' plots show the difference between nt(r) and the 
fit with n t (r) fit . 



sible task since the result is unavoidably affected by the 
precision on the sphere centres and the polydispersity of 
the spheres themselves. With the data from X-ray to- 
mography, we can calculate the location of the sphere 
centers with a precision which is well within 1% of their 
diameters. On the other hand, the beads utilized have a 
polydisperse bead-diameter distribution with a standard 
deviation around ~ 0.02d Therefore, statistically, the 
large majority of neighbors in contact must stay within a 
radial distance of 1.02d. Table[I]reports the values of the 
average number of neighbors (n t ) computed in G at the 
four different radial distances: r — d, 1.02<i, 1.05d and 



FIG. 4: In the region r < d the distribution of radial dis- 
tances among neighboring spheres are very well mimicked by 
a Normal distribution (Ant(r) being the average number of 
spheres at radial distance between r — Ar and r + Ar) . All 
data for all the 6 samples collapse into a single behavior when 
Ant(r)t; is plotted versus (r/d— l)/£. The line is the (col- 
lapsed) fit of the data in the region r < d (left-hand side) 
with the Normal distribution. This is the same fit used in the 
deconvolution shown in FigE] which allows to estimate the 
average number of spheres in contact. 

l.ld for the 6 samples A-F. We observe values for the av- 
erage number of neighbors between n t (r) ~ 2.9 and 8.9, 
and an increasing trend with the packing density. 



B. Touching Neighbors 

A more precise estimate for the actual number of 
spheres in contact can be inferred from the behavior of 
n t (r) (shown in Fig|5] as function of the radial distance 
up to r = l.ld). From FigO one can note that above 
r = 0.98d the number of neighbors grows very steeply up 
to a 'knee' at about 1.02d where a slower growth takes 
place. Such a steep growth in the number of neighbors 
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FIG. 5: Number of neighbors in contact vs. sample density. 
The filled symbols corresponds to the samples investigated in 
the present work. The two symbols ' + ' are the values from Q 
whereas the ' x ' is from [2, ■ The line is the best-fitting of 
the A-F data with a linear trend constrained to pass through 
n c = 4 at p — 0.55. 



can only be an effect of the uncertainty on the positions 
of the sphere centers and of the spread in the statistical 
distribution of the distances between spheres in contact 
(which is a consequence of the polydispersity). Indeed, 
in the ideal case when all the exact positions of closely 
packed, identical, perfect, spheres are known, one would 
expects that nt(r) has a discontinuity at r — d (from 
zero to the number of neighbors in contact n c ) followed 
by some kind of growth for r > d. For real, polydisperse, 
non-perfect, spheroidal grains, the distance between el- 
ements in contact is not a fixed value but instead it is 
distributed around an average value. In Fig 0] it is clearly 
shown that such a distribution of radial distances is well 
mimicked for r < d by a normal distribution. As a conse- 
quence, we expect that the steep growth of nt{r) around 
r = d is well described with a complementary error func- 
tion (for r < d). To such error function we must add the 
contribution from the 'nearly touching' spheres which is 
expected to become sizable from r > d. We therefore 
expect to find an error-function behavior up to r ~ d 
and then a combined contribution from the error-function 
and some growing law describing the cumulate number 
of non-touching neighbors within the distance r. Indeed, 
we verify that the behavior of rit (r) for r < d is very well 
described by a complementary error function normalized 
to n r : 



n t (r) 



fn 



1 



eX P( ^2~> dx 

i Z S 



(1) 



Where the value of the mean (d) is the average sphere 
diameter which was estimated: d — 25.00 voxels (samples 
B D E F) and d = 30.81 voxels (samples A C) 0. On the 
other hand, the variance £ can be directly estimated from 
the data by computing the second moments ((r — d) 2 } for 
the radial distances between spheres calculated over half 
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with i,j indices labeling the sphere centers; the symbol 
rij indicates the distances between the centers of sphere 
i and j; and H(d— rij) the step function which returns 
1 if Tij < d and if r^j > d. From Eq[21we retrieve vari- 
ances £ in the interval between O.Old and 0.02d (all the 
values are reported in Table [I] ) . These values are con- 
sistent with the bead-polydispersity and are significantly 
larger than the estimated uncertainty on the sphere cen- 
ters. The only free parameter left in Eq^ is the value 
of n c which can be now computed by best fitting the 
agreement between the data for n t (r) and the function 
n£ l *(r) in the region r < d. In Fig[3|it is shown that the 
function nl lt (r) fits well the data for r < d by using the 
values of n c given in Table [I] At larger distances (r > d) 
near-neighbors not in contact start to contribute to nt (r) 
and the two behaviors split. We estimate that in the 6 
samples A-F there are in average between 5.81 and 6.97 
spheres in contact (see Tableland FigEJ). These num- 
bers fall in the range of reported values: Bernal measures 
n c = 5.5 at p = 0.6 and n c = 6.4 at p — 0.62 p); whereas 
from the data by Scott we have n c = 7 at p = 0.63 
. Note that this latter data has been recalculated from 
that reported by applying the deconvolution method de- 
scribed above. 

In Fig [5] the values of n c vs. density for the samples 
A-F are reported together with those from Bernal and 
Scott. Such agreement between these different data is 
remarkable considering the different experimental proto- 
cols, the different preparations of the samples, the differ- 
ent criteria for identifying and counting spheres in con- 
tact and the different polydispersity of the spheres. As 
one can see from Fig03 they all show a clear and con- 
sistent increasing behaviour with the density. A similar 
increasing trend was also found in simulated packings 
|20l 12 1| . This dependence on the packing density has 
important theoretical implications which are discussed 
hereafter. 



Mechanical Equilibrium 

In a stack of grains at mechanical equilibrium, New- 
ton's equations for the balance of force and torque acting 
on each grain must be satisfied. Lagrange and Maxwell 
[22I |23| have been the first to note that in these kinds 
of systems, to achieve stability, the number of degrees 
of freedom must balance the number of constraints. It 
is straightforward to calculate that the balance between 
freedom and constraints requires n c = 6 in the case of 
perfectly spherical friction-less spheres, and n c = 4 for 
more realistic grains (non-spherical with friction) 
These values of n c are encouraging. Indeed, they predict 
values which are rather close to the ones observed exper- 
imentally. However, it must be noted that this condition 
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for n c is neither sufficient nor necessary |26| |27J, [2j| . In- 
deed, there can be local configurations which contribute 
to n c but do not contribute to the rigidity of the whole 
system. (These are, for instance, the 'rattling' grains 
which can be removed from the system without affecting 
its stability.) On the other hand, there are local arrange- 
ments which satisfy the counting rule on n c but never- 
theless are not rigid 26, 27]. 

In recent years there have been scores of theoretical ap- 
proaches which consider real, disordered, granular pack- 
ings to be isostatic (free of self induced stresses) [2J, |22j . 
The advantage is that in a system at isostatic equilib- 
rium, the intergranular forces are uniquely determined by 
the balance of force and torque alone. On the contrary, 
a overconstrained structure can generate self-stress and 
the deformation of individual grains becomes relevant. 
In real granular materials (or in bead packs) friction and 
rotational degrees of freedom are unavoidable, therefore 
the Maxwell counting implies that isostatic configura- 
tions must have average connectivity of n c = 4. Unfor- 
tunately, such a value is rather small in comparison with 
all available empirical estimations. Moreover, all the ex- 
perimental observations conclude that n c increases with 
the packing density excluding therefore the possibility to 
fixing n c at 4 for all sTable packings. 

However, in related studies for the rigidity in network 
glasses [2^, Hcj it has been observed that there can ex- 
ist two phase transitions associated with the increase of 
connectivity in a network: a rigidity percolation and a 
stress percolation, and between these two thresholds an 
intermediate phase which is rigid but unstressed. The 
rigidity percolation occurs at the threshold predicted by 
the Maxwell counting (n c = 4) and appears to be a sec- 
ond order of transition |3C| . The stress percolation tran- 
sition occurs at an higher value of n c and it could be 
a first-order transition |3(J. This suggests that granular 
packings might be in a marginal state between the rigid- 
ity and the stress percolation thresholds. In such an iso- 
static unstressed state, the system has zero elastic mod- 
ulus (in the thermodynamic limit) [2^, [3(| , it is therefore 
marginally rigid and it can be seen as in a intermediate 
state between fluid and solid ■ An extrapolation from 
the experimental data for n c , reported in Fig[5l suggests 
that a grain connectivity equal to 4 could be reached by 
the system at the density p — 0.55. This would place the 
rigidity percolation threshold at the loose packing limit 




FIG. 6: Same as Fig^with a portion removed and the topo- 
logical distances from a given central sphere highlighted in 
colours (online version). 
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FIG. 7: Shell occupation numbers vs. topological distance. 
The symbols indicate the different samples (as in Fig|3J and 
the lines are the best-fits (of the growing part only) using the 
polynomial form: Kj = aj 2 + aj + Co. The fits are between 
3 = 2 and 3 = 10 (for samples B, D, E, F ) and j — 15 (for 
samples A, C). The data refer to threshold distances 1.02d, 
1.05d, Lid and lAd as indicated in the figures. 



IV. CONTACT NETWORK: BEYOND FIRST 
NEIGHBORS 

Any force-path or any infinitesimal local grain dis- 
placement must mechanically propagate from grain to 
grain through the network of touching grains. The un- 
derstanding of the hierarchical organization of such con- 
tact network beyond first neighbors is therefore of great 
importance. Here we apply to granular matter an ap- 



proach which was orig inally developed for the study of 
crystalline systems |33l 134 l35l l36t I37L l38| and disordered 
foams 0, 113, |4jj- The topological structure of crys- 
talline frameworks has been intensely studied in terms 
of the number of atoms that are 7 bonds away from a 
given atom [H H HI IM EUSi . If we start from a 
given 'central' atom, the first 'shell' (distance j = 1) is 
made by all the atoms in contact with the central one. 
The second shell (distance j = 2) consists of all atoms 
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FIG. 8: The coefficient a increases with the threshold 8 on 
the radial distance (<5 is expressed in sphere-diameters unit). 
The insert shows that, for all the samples, the coefficient a 
follows the linear law a ~ b(r — d) + ao in the region between 
1.05 and 1.1. 

which are neighbors to the atoms in the first shell, ex- 
cluding the central one. Moving outward, the atoms at 
shell 3 + 1 are all the ones which are bonded to atoms in 
shell j and which have not been counted previously. In 
infinite, periodic, crystalline structures with no bound- 
aries, the number of atoms per shell should increase with 
the topological distance and it has been shown that in 
several three-dimensional crystalline structures the law 
of growth for the number of atoms (Kj) at shell j can be 
described with: Kj = ajj 2 + bjj + Cj (with a,j, bj and 
Cj coefficients that might vary with j but only within 
a bounded finite interval) Hi 1^1 El IH S E3- 
Following the definition of O'Keeffe 36], for these crys- 
talline systems, the asymptotic behavior of Kj can be 
characterized in terms of an 'exact topological density' : 
TD = (a,*) /3 [41), It has been noted that such a topo- 
logical density is interestingly related to the geometrical 
density of the corresponding crystalline structure and it 
is a powerful instrument to characterize such systems. 
For instance, it is easy to compute that the cubic lat- 
tice has Kj — 4j 2 + 2. Whereas, spheres packed in a 
bec (body centered cubic) crystalline arrangement have: 
Kj = 6.7 2 + 2 (j > 0). On the other hand, it has been 
shown [38j that for Barlow packings of spheres, Kj are 
always in a narrow range within: 

21 i 2 

10f +2<Kj< +2 (j > 0) , (3) 

where the brackets [—J indicate the floor function. In 
EqOHthe lower limit is associated with the fee (face cen- 
tered cubic) packing and the upper limit corresponds 
to the and hep (hexagonal closed packed) packing. It 
has been observed by O'Keeffe and Hyde jU that for 
lattice sphere packings the general rule holds: Kj — 
(n c — 2)j 2 + 2, implying therefore a — n c — 2. 

Beyond perfect crystalline order very few results are 
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FIG. 9: The coefficient a plotted against the average number 
of neighbors in contact n c show that disordered packings have 
larger topological densities in comparison with lattice sphere 
packings. 



known either from theoretical, empirical or numerical 
point of view. One can argue that Kj must grow with 
a law comparable with the law for a spherical shell: 
Kj ~ aj 2 ~ 4-7r/ 2 . However, it is also clear that the 
shape of the growing shell and its roughness can dras- 
tically change the coefficient a (as observed in two di- 
mensional cases [3!j EH)- Moreover, it can be shown 
[44| that in some topological networks the law of growth 
can follow an intrinsic dimension which is different from 
the dimension of the embedding space (3 in our case). 
This mechanism can produce power law growth with ex- 
ponents different from 2, or different behaviors such as 
exponential - or even faster - laws of growth |44| . 

We observed that the number of spheres at a given 
topological distance j from a central one follows a power 
law growth (see Fig0 until a critical distance j, above 
which the shells hit the sample boundaries and Kj starts 
to decrease. We verify that a quadratic law Kj = 
a J 2 + c ii + c o fits quite accurately the observed behav- 
iors of Kj for j < j. This fixes the intrinsic dimension 
for these systems equal to 3 (which coincides with the 
geometrical dimension of the embedding space). The co- 
efficient a depends on the threshold 5 on the radial dis- 
tances within which we consider spheres to be connected 
in the contact network. Indeed, changes in the thresh- 
old distances are unavoidably associated with changes in 
the contact network and an enlargement of the threshold 
distance must correspond to a thickening of the shell. In 
FigJSl we show that the dependence of the coefficient a 
on the threshold distance 6 is rather complex and not 
reducible to a simple law. However, we verify that in the 
interval of threshold distances between 1.05d and l.ld, a 
simple linear increment is observed: a ~ bd + ao (typi- 
cally with b ~ 0(1)). This is shown in the insert in Fig|8] 
Such a law suggests that a unique value for the topologi- 
cal density, independent from the threshold distance, can 
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be associated with the value of the coefficient a at 5 = 0: 
T D = a /3. In Table UTI are reported the topological den- 
sities for all the samples A-F. As one can see the quantity 
3Td = do stays in a rather narrow range around 8.5 and 
slightly increases with the sample density. Interestingly 
we observe in FigElthat disordered sphere packings have 
coefficient a consistently above n c — 2 implying therefore 
that such packings have larger topological densities than 
lattice sphere packings with the same coordination num- 
ber. This observation might be relevant when the struc- 
tural stability and rigidity of such system is concerned. A 
view of the topological shell structure constructed from 
a given central sphere in one of the samples, is shown in 
Fig© 



V. LOCAL ORIENTATION 

Revealing and quantifying orientational order is a key 
issue in establishing the nature of internal organization, 
and in particular in determining whether there exists a 
'typical' disordered state or identifying possible tenden- 
cies towards hidden symmetries. Indeed, if such a 'typ- 
ical' state exists or/and if there is a tendency toward 
a specific local organization, then it will be possible to 
associate to a given granular pack an order parameter 
which could measure how close the packing is to the ideal 
structure. On the other hand if one can prove that the 
system is a collection of uncorrelated local configurations 
then this will make it possible to calculate the configura- 
tional entropy and -consequently- the probability to find 
the system in a given state at a given density. It has 
been often argued that the competition between the ten- 
dency to form a locally compact configuration and the 
geometrical frustration could be the key to understand 
the mechanism of formation of disordered packings and 
glassy structures. If this is the case we will expect to 
see at local level, configurations with rotational symme- 
tries characteristic of icosahedral and other closed packed 
structures. The study of the local rotational symmetry 
can therefore give insights also on the mechanism of for- 
mation of these structures. 

The challenge is to find a measure of rotational sym- 
metry which is invariant with respect to rotations in the 
system of coordinates. A powerful solution was intro- 
duced by Steinhardt, Nelson and Ronchetti by as- 
signing a set of spherical harmonics Yj, m (0(r), (f>(r)) to 
the vectors r between couples of spheres (with 9(r) and 
4>{r) the polar and azimuthal angles of r) and introducing 
the quantities: 



A-7T 

\ wi=—l 



(4) 



with average ((...)) over the bonds in the sample. 
Such a quantity is invariant under rotations in the co- 
ordinate system and it takes characteristic values which 
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FIG. 10: Values of (Qi,Q6) for all the local configurations 
in samples 'A-F' given a threshold distance of l.ld. Each dot 
correspond to a given sphere in G. The color (online version) 
ranges from light-red to dark-red depending on the number of 
neighbors (within the threshold distance) of each local config- 
uration. The lines are contour plots of the frequencies. The 
positions of specific symmetries (ico, sc, bcc, fee, hep) in the 
(QiiQe) plane are also indicated with '*' and projected on 
the axes with dotted lines . 



can be used to quantify the kind and the degree of ro- 
tational symmetry in the system. However, it must 
be noted that the quantity Yi m (6{?i), 4>(fi)) depends 
on the orientation, therefore in the case of a polycrys- 
talline aggregate, with finite correlation length, its av- 
erage (Yi ym (9(fi),(j)(fi))} will decrease and tend to zero 
with the sample-size. To avoid this inconvenience, which 
makes the comparison between values of Qi on differently 
sized samples meaningless, it is convenient to adopt a lo- 
cal measure of Qi by restricting the average only over 
the local bonds between a sphere and its neighbors. In 
this way, to each sphere in the system can be associated 
a Qi and local order can be singled out by counting the 
number of configurations with Qi corresponding to spe- 
cial symmetries. In particular the two cases I = 4 and 
I = 6 have special significance. For instance, the simple 
cubic lattice has (Qi,Q e ) sc = (0.764,0.354), the body 
centered cubic lattice has (Q4,Qe) bcc = (0.036,0.511), 
the fee has (Q 4 ,Q 6 ) /CC = (0.191,0.574), the hep has 



\hcp 



(0.097, 0.485) and the icosahedral rota- 



tional symmetry gives {Qi, Qe) lco = (0, 0.663). Since the 
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3T D 


r treshold 


(Qi,Qe) 


dis (%) / 


sc (%) hep (%) 




1.1 


(0.27,0.47) 


23 


3 


1 


A 7.2 ±0.3 


1.2 


(0.22,0.42) 


32 


2 


3 




1.3 


(0.18,0.40) 


38 


1 


5 




1.4 


(0.15,0.36) 


42 
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4 




1.1 


(0.30,0.45) 


24 


3 
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B 7.2 ±0.4 


1.2 


(0.23,0.44) 


32 
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3 




1.3 


(0.16,0.38) 


37 
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5 




1.4 


(0.14,0.35) 


43 
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5 




1.1 


(0.23,0.46) 


28 
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C 8.7 ±0.4 


1.2 


(0.21,0.43) 


35 
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7 




1.3 


(0.15,0.40) 


41 
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11 




1.4 


(0.12,0.37) 


45 
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1.1 


(0.25,0.44) 


28 
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D 8.4 ± 0.3 


1.2 


(0.19,0.44) 


35 
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7 




1.3 


(0.15,0.40) 


42 
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11 




1.4 


(0.11,0.36) 


46 


1 
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1.1 


(0.22,0.44) 


27 


5 
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E 8.6 ±0.4 


1.2 


(0.20,0.43) 


37 
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7 




1.3 


(0.15,0.39) 


42 
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12 




1.4 


(0.12,0.36)) 


47 
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1.1 


(0.23,0.44) 


31 


6 


4 


F 8.9 ±0.4 


1.2 


(0.16,0.45) 


38 
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12 




1.3 


(0.13,0.42) 


43 
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17 




1.4 


(0.10,0.38) 


47 


3 
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TABLE II: Topological Densities (3Td = ao). Most recurrent values for the local orientation order (Qi, Qe)- Fraction of local 
configurations with (Qi, Q§) in the range (Qi ± 0.05, Qe ± 0.05) (dis). Fraction of local configurations close to special form of 
order: (fee) with (Qi,Q 6 ) in the range (0.191±0.05, 0.574±0.05) ; (hep) with (Qi,Q 6 ) in the range (0.097±0.05, 0.485±0.05) 



lowest non-zero Qi common to the icosahedral, hexagonal 
and the cubic symmetries is for I = 6, it has been argued 
by several authors that the value of Q 6 is a good indica- 
tor of the degree of order in the system and it might be 
used as an 'order parameter' 0, EE EH EE EH Indeed, 
Qe is very sensitive to any kind of crystallization and it 
increases significantly when order appears |45| . 

Similarly to that discussed in the previous sections, the 
measure depends on the adopted geometrical criteria to 
identify neighbors. In the literature, several different cri- 
teria are used: in ^2] ai l neighbors within 1.2d are consid- 
ered; in the neighbors up to the radial distance which 
correspond to the first minimum in the radial distribu- 
tion functi on ( r ~ lAd, see section IYlI|) was considered; 
in [45| and [43 the Voronoi (or Delaunay [lj, |5(j ) neigh- 
bors where used instead. This last definition might be 
misleading in some cases (as pointed out by |l2j|). since 
the Voronoi method tends to associate bonds to some- 
time distant neighbors. For instance, an fee crystalline 
arrangement (with infinitesimal perturbation) takes 2 ex- 
tra neighbors (from 12 to 14 in average) using the Voronoi 
criteria. Here the influence of the neighboring criteria is 
analyzed by using 4 different threshold distances: l.ld, 



1.2d, 1.3d and I Ad. 

Examples of the distribution of local (Qa,Qq) are 
shown in FigsllUlandllll 

We observe values of (Qa,Qq) narrowly distributed 
around their most recurrent values (Qa,Q§) with very 
large fractions of local configurations (between 23 and 
47 %) which have local symmetries characterized by 
(Q4,Qe) within the range (Q 4 ± 0.05, Q 6 ± 0.05). The 
values of (Q&,Qs) range between 0.10 < Q4 < 0.30 and 
0.35 < Qe < 0.45 across all samples and all thresholds 
(see Table ITlTl . Such values are far from any special sym- 
metry. In order to search for signatures of known local 
symmetries we measured the fraction of local configu- 
rations with (Qi, Qe) in a region within a range ±0.05 
from the values in the ideal structures ( fee, hep, ico, sc, 
and bee). We found that there are no significant frac- 
tions (below 1 %) of local configurations with symmetry 
compatible with icosahedral, simple cubic or bee; there 
is a small fraction (between 1 to 6 %) of configurations 
with local symmetry compatible with fee, and there is 
a fraction of configurations with /icp-kind of local order 
which becomes quite significant at large densities (reach- 
ing 17% at p = 0.64 and S — 1.3d). This occurrence of a 
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FIG. 11: Same plot as in Fig llOl but with threshold distance 
lAd. 



rather large fraction of local symmetry with an hcp-\ike 
character might suggest the beginning of a crystallization 
process. However, we have verified that there are no cor- 
relations between neighboring sites with symmetry close 
to hep. This excludes the presence of any long range hep 
order recurrent or symmetrical organizations beyond first 
neighbors. 

These findings cast considerable doubts over the exis- 
tence of any crystalline order and also question the idea 
that a tendency toward local -frustrated- icosahedral or- 
der can be responsible for the resilience to crystallize for 
such packings. These results will be confirmed and rein- 
forced in Section IVlIII where the local densities are stud- 
ied. On the other hand, these findings are not conclusive 
because the origin and nature of the most abundant con- 
figurations with {Qi,Qe) ~ (0.25,0.45) are still elusive. 
Further studies to clarify the nature and the origin of 
such local configurations are needed. 



VI. PACKING EFFICIENCY 

In this section we investigate how the global sample 
density is perceived by a sphere at the local level and how 
the global packing affects the local environment. To this 
purpose we compute the number of spheres placed within 
a certain radial distance from a given sphere. This quan- 
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FIG. 12: Total number of sphere centers within a given ra- 
dial distance r: average, maximum and minimum for all the 
samples. The filled area is bounded at the top by the most 
efficient packings among the fee, hep, bee, sc. 



tity (n t (r)) can be viewed as a measure of how efficiently 
local dense agglomerates of spheres are formed. We called 
this measure the packing efficiency 0. It is well known 
that no more than 12 spheres can be found in contact 
with one sphere (the 'Kissing number' [l3lp. but the up- 
per limit for the number of spheres within a given radial 
distance is, in general, unknown. Fig ll2l shows the aver- 
age, the maximum and the minimum numbers of neigh- 
bors within a given radial distance from any sphere in G. 
Clearly, there are no neighbors up to distances close to 
r ~ d, when suddenly the number of neighbors increases 
very steeply and then, after this jump, it increases with 
distance following a less steep trend with very compara- 
ble behaviors between all 6 samples A-F. It is of some 
interest to compare the values of nt(r) empirically ob- 
tained for these disordered samples with the known ones 
associated with crystalline structures. We observe that 
in a large range of radial distances between Id and 2d 
there are some local configurations with packing efficien- 
cies which are above the crystalline ones. Moreover, in 
the region around r ~ 1.6d, disordered packings show 
better average packing efficiencies than the crystalline 
ones (fee, hep, bee, sc). This is rather surprising if we 
consider that the fee and hep packings are more than 
15% denser than the disordered ones. However, one can 
note that increasing the distance, the disordered pack- 
ings become less and less efficient in comparison with 
the close packed crystalline arrangements, and above 2d 
all the configurations have a smaller cumulate number 
of neighbors than the close crystalline packings. Inter- 
estingly, in the same region around 1.6d where the dis- 
ordered packings are very efficient we also observe the 
minimum spread for the values of the efficiencies across 
the samples at different densities. We discuss in further 
details the behavior of the number of neighbors with the 
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FIG. 13: (a) Normalized radial distribution function, (b) The 
detail of the two peaks respectively at y/2d and 2d (vertical 
lines), (c) The ratio between the value of the peak at 2d and 
the one at 



radial distance in the next Section where the radial dis- 
tribution function is analyzed. 



VII. 



RADIAL DISTRIBUTION FUNCTION 



The Radial Distribution Function (g{r)) is the prob- 
ability distribution to find the centre of a particle in a 
given position at distance r from a reference one. This 
measurement is widely used in geometrical characteriza- 
tion of packing structures and contains information about 
long range interparticle correlations and their organiza- 
tion UllS 

In order to calculate this quantity one must count the 
number of sphere centers within a radial distance r from 
a given sphere centre. The average of this number is the 
quantity nt(r) studied in the previous section and it is 
related to the radial distribution function by: 



Thin) - n t (r ) 



g(r)4:TTr 2 dr 



(5) 



Therefore, given the position of the sphere-centers, these 
two quantities rit(r) and g(r) can be straightforwardly 




FIG. 14: The three peaks of the Radial distribution function 
respectively (from top to bottom) at r = 1, \/3 and 2 can be 
described with power law singularities: g(r) ~ co\r — ro\~ a . 
The coefficient a depends on the sample-densities and their 
behaviors are reported in the inserts. 



computed. Here, we calculate the normalized radial dis- 
tribution function g(r) which is the average number of 
sphere centres, within a radial distance r — A/2 and 
r + A/2, divided by cr 2 . With the constant c fixed by 
imposing that asymptotically g(r) — > 1 for r — > oo. (We 
have verified that different choices of A within a broad 
range of 10 _4 d to 10 _2 <i lead to almost indistinguishable 
results.) 

In Fig. the behavior of g(r) vs. r/d is shown. We 
observe a very pronounced peak at r = d which corre- 
sponds to the neighbors in contact. Then the probability 
to find neighbors decreases with r reaching a minimum 
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around lAd. Subsequently, at larger radial distances, 
the probability increases again forming two peaks respec- 
tively at r = V3d and r ~ 2c? and then after these peaks 
it continues to fluctuate with decreasing amplitudes. The 
details of the second and third peaks, plotted in Figll3b. 
show that the two peaks at r = \/3d and r ~ 2d both in- 
crease in height with the packing density. With the peak 
at r — 2d growing faster than the one at r = V3d (see 
Figll3bh This might indicate an increasing organisation 
in the packing structure but, on the other hand, no signs 
of crystallisation were detected (see, Section Section 

EnH). 

For all the samples investigated, we found that the 
behavior of g(r) at radial distances between r ~ Id and 
r ~ lAd (between the first peak and the first minimum) 
can be quite accurately described in terms of a power law 
singularity: 

g( r ) ~ : C " : , (6) 

V 1 \r-r \ a y ' 

with good fits for ro = 1.03 and a which increase with 
the sample-density from 0.27 to 0.45 (Fig |14fa h A similar 
behavior, but with a — 0.5 and ro = d, was reported in 
|19| for numerical simulations. A more recent numeri- 
cal investigation proposes an exponent a ~ 0.4 [52^. In 
Fig ll4l it is also highlighted the growing trend of a with 
the density p. 

Interestingly, the behavior of g(r) around the other two 
following peaks (respectively at r ~ y3d and ~ 2d) can 
be described by using similar power law kind of diver- 
gences. In particular the region 1.4 < r < 1.73 is well 
fitted by Eq|fj|with ro = 1.8 and a between 0.37 and 0.62. 
Whereas, the region 2 < r < 2.15 is well fitted by using 
ro = 2 and a between 0.11 and 0.17 ('Fig ll4|l . We must 
stress that these are qualitative behaviors: a reliable fit 
with a power law trend must be performed over several 
orders of magnitude in the x and y scales. These linear 
interpolations in log-log scales must therefore be consid- 
ered more as indicative behaviors of qualitative laws more 
than fits. 

The origin and the nature of this power-law like be- 
havior around these peaks is rather puzzling. Indeed al- 
though the presence of such peaks clearly indicate some 
kind of organization in the system, on the other hand 
other analyses, such as the orientational symmetry dis- 
cussed in Section exclude the presence of any crys- 
talline or polycrystalline pattern in the samples. To bet- 
ter understand this issue one must single out the specific 
organization of the local configurations which contribute 
to each peak. This will be the topic of a future paper. 

VIII. GLOBAL AND LOCAL DENSITIES 

We have already referred previously to the sample den- 
sity which is the fraction between the volume occupied 
by the spheres divided by the total volume. Although, 
this definition is very straightfoward, on the other hand 
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FIG. 15: (above) Distribution of the local densities in G 
for the 6 sapmles. (below) The cumulants show tails that 
decrease exponentially or faster with slightly asymmetric dis- 
tributions in the left and right parts. 



it is well defined only for an infinite sample. In all the 
other practical cases, where boundaries are present, the 
density is unavoidably associated to the way of partition- 
ing space. A convenient way to study the density is by 
partitioning the space in local portions and introducing a 
local- density associated with the fraction of volume occu- 
pied within each local portion of space. Surprisingly, in 
the literature of granular matter, very few investigations 
have been devoted to the study of local densities either in 
experiments or in simulations. On the contrary the un- 
derstanding of how the space is shared among the packed 
spheres and finding how efficiently the spheres can pack 
locally is an essential information which can contribute 
to the understanding of both the structure and the for- 
mation of these systems. 

We calculate the local- densities which are defined as 
the fractions between the sphere-volumes and the vol- 
umes of the Voronoi cells 50] constructed around the cen- 
tre of each sphere in the sample (recall that the Voronoi 
cell is the portion of space closest to a given centre in re- 
spect of any other centre) . The sample-densities are frac- 
tions between the sum over the volumes of the spheres 
in G and the sum over the volumes of the Voronoi cells 
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associated with these spheres. We observe that typically 
the density is not homogeneously distributed in different 
parts of the samples. This is also discussed by several 
other works (sec for instance jHJ I n the samples 

A-F the densities are relatively smaller than the average 
in a region close to the cylinder central axis; the density 
increases moving outwards from the centre, then it sat- 
urates to rather homogeneous values up to a distance of 
a few (2-3) sphere diameters from the boundary. Rather 
inhomogeneous densities are also observed in the vertical 
direction, but in this case we find different behaviours 
depending on the sample-preparation. However, we ver- 
ify that in all the samples sub-regions the densities stay 
in a rather narrow range (within ± 0.01, see Table HJ 
form the average ones. More importantly we verify that 
all the computed structural properties do not change sig- 
nificantly in their behaviors and characteristics with the 
part of sample analyzed. In the second column of Ta- 
ble |U the average density values for samples A-F and 
their interval of variations are reported. 



Local Density distributions and geometrical 
frustration 
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FIG. 16: Probability to find a path through the first neigh- 
bors with section larger than g. 

less than 0.07% for A and B (p = 0.586 and 0.593) This is 
a very strong indication that local sphere-arrangements 
with high local densities -such as the icosahedral packing 
- play no role in these disordered sphere packings. 



Figure IT51 shows the distribution of the local densities 
in G for the 6 samples. We observe that these local den- 
sities have slightly asymmetric distributions with expo- 
nential -or faster decays from the average densities (which 
are in the range 0.586 < p < 0.640) and have standard 
deviations a within 1.5 %. 

It has been often argued that the driving mechanism 
which generates amorphous structures could be the geo- 
metrical frustration. This derives from the fact that the 
locally densest configuration in equal sphere packing is 
achieved by placing on the vertices of an icosahedron 12 
spheres in touch with a central sphere. But, iocosahe- 
dral symmetry is incompatible with translational sym- 
metry and this generates frustration: some gaps must be 
formed and the symmetry must be broken. Indeed, it is 
known 0] that although there exist locally denser con- 
figurations, at global scale the densest achievable pack- 
ings have density p < 7r/\/l8 = 0.74048... which is the 
one of fee and hep crystalline packings. 

Here we test whether such a frustration mechanism 
is really relevant in our granular packs. To this end 
we search for local configurations which are locally close 
packed at densities larger than the crystalline ones. If a 
sizable amount of such configurations is found this im- 
plies that indeed the geometrical frustration must have 
an important role in amorphous systems. The result is 
unexpected: only 14 local configurations with densities 
above 0.7405 were found over a set of more than 209,000 
local configurations in the 6 different samples. Moreover, 
very few local configurations have density between 0.7 
and 0.7405. Respectively we find: less than 1.7 % in the 
densest sample (F at p = 0.640); less than 1% in the three 
samples C, D, E (with p = 0.617, 0.630 and 0.626); and 



IX. GEOMETRY AND STRUCTURAL ARREST 

We calculated the gaps between neighboring spheres 
and evaluated the size of the largest gap for each local 
configuration surrounding a given sphere. This gives the 
probability for each sphere to move outward from a given 
local configuration without displacing the positions of its 
neighbors. Such an 'escape probability' is computed by 
constructing circles through the centres of the Delaunay 
simplexes incident on the central sphere. The central 
sphere can 'escape' from this configuration only if its di- 
ameter (d) is smaller or equal to the largest radius of such 
circles. Clearly, when such a move is possible, the system 
can change its geometrical configuration by means of lo- 
cal moves only. From a thermodynamical point of view 
this implies that it can dynamically explore the phase- 
space with low energy moves and reach equilibrium in 
short relaxation times. On the contrary, when the es- 
cape probability is zero, a structural rearrangement re- 
quires the displacement of a larger number of spheres and 
the system is more likely to be trapped for long times in 
meta-stable states. 

The relative number of local configurations with gaps 
larger or equal to a threshold size g is reported in Fig ll6l 
for the 6 samples A-F. We find that all the samples with 
p > 0.6 (C-F) do not have any configuration which allows 
the central sphere to 'escape'. Whereas samples A and 
B have few local configurations with gaps larger than d 
but seems statistically irrelevant (< 0.1%). This strongly 
suggests that around p ~ 0.58—0.6 an important phase in 
the system dynamics reaches an end: above this density, 
local readjustments involving only the displacement of a 
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single sphere are forbidden. The particle mobility is con- 
strained mostly within the Voronoi cell and the system 
compaction can proceed only by involving the collective 
and correlated readjustment of a larger set of spheres. 
At this stage the system can no longer sample the whole 
phase-space and it is trapped within the basin of attrac- 
tion of some inherent configuration |55| and eventually 
will reach a structural arrest before the thermodynami- 
cal equilibrium is reached. 

X. CONCLUSIONS 

The structure of disordered packing of mono-sized 
spheres has been investigated by means of X-ray Com- 
puted Tomography. We performed an extensive study 
over 6 large samples at packing densities ranging be- 
tween 0.586 and 0.640, investigating several geometrical 
and topological properties. 

The number of neighbors surrounding each sphere in 
the packing was studied with unprecedented statistical 
accuracy fScction lIII|) . The average number of spheres in 
contact was extracted by means of an innovative method 
which deconvolutes the contribution of touching neigh- 
bors from the contribution from near neighbors. The re- 
sults show that the average number of spheres in contact 
increases with the sample density and it is between 5.5 
and 7.5 in the range of densities examined (see FigJSJ. 
An extrapolation to the random loose packing density 
(p = 0.55) suggests that at this density the system could 
have an average number of 4 neighbors per sphere. This 
implies the possibility of a rigidity percolation transition 
taking place at the random loose packing limit. 

The structure beyond first neighbors, studied by means 
of a topological map ( Section Hv|) . shows that the con- 
tact networks of these systems have intrinsic dimension 
3 and topological densities between 7.2 and 8.9 (see Ta- 
ble [HJ . This is a novel approach for three-dimensional 
non-cystalline structures. Surprizingly, we found that 
the topological density in disordered sphere packings is 
always larger than the topological density in the corre- 
sponding lattice sphere packings. Such a larger topo- 
logical density is an indication that the contact network 
is more compact in disordered systems despite the fact 
that the geometrical density is lower. This fact might 
have important implication on the system stability and 
resilience to perturbations and shocks. 

A non-intuitive property was found by computing the 



packing efficency ( Section I Vlfl . It results that disordered 
packings can have a larger number of neighbors within 
a given radial distance than the crystalline packings (see 
Fig ll2f> . This fact is surprising if one considers that fee 
and hep packings are more than 15% denser than the 
disordered packings. 

In Section [V] we searched for local symmetries and 
the existence of significant repetitive local configurations. 
The result was that most of the spheres are arranged 
locally in configurations which are significantly different 
from any crystalline arrangement excluding therefore the 
presence of any partial crystallization in these samples. 
We also established that there are no statistically signif- 
icant configurations with icosahedral symmetry. More- 
over in Section lVIIII wc analyze the local densities of more 
than 200,000 configurations concluding that there are no 
statistically significant local arrangements with density 
equal or above 0.74. This excludes that geometrical frus- 
tration can play any significant role in the formation of 
such amorphous packings. 

The structural organization emerging from the radial 
distribution function was discussed in Section I V 111 where 
we pointed out peculiar power law kind of behavior 
around the peaks of g(r) at r = d, r = ^/2d and r = 2d 
(see EqEJand FigHU). 

Insights on the dynamical formation of these systems 
were given in Section llXI where we pointed out that the 
distribution of gaps between neighboring spheres sug- 
gests that a dynamical glass transition might take place 
at densities around 0.58-0.60. 

With the largely increased statistical confidence, these 
datasets may open new and more fruitful paths to the 
understanding of granular materials. 
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